#ifndef _VEC_H_
#define _VEC_H_

#include<iostream>
#include <utility>
#include <cmath>
#include <cassert>
using namespace std;
// todo : support for constexpr

/// Define a point or a vector in the Euclidean space.
/**
 */
template <class T1, int _Dim>
class Vec
{
protected:
  T1 coord[_Dim];

public:
  enum { Dim = _Dim };
  Vec(){};
  Vec(T1 a[]){
  for(int i=0;i<Dim;i++)
	coord[i]=a[i];
   }
  
 

  T1 &operator[](int _d) { return coord[_d]; }

  const T1 &operator[](int _d) const { return coord[_d];  }

  const T1 *data() const { return &coord[0]; }
  void set(T1 x,int i){coord[i]=x;}




template <class T>
inline Vec<T,2> clockwise(const Vec<T,2> &v) {
  return Vec<T,2> {v[1], -v[0]};
}

//==================================================
// operations on multi-dim vectors

template <class T, int Dim>
inline Vec<T,Dim> operator-(const Vec<T,Dim> &lhs) {
  Vec<T,Dim> res;
  for(int d=0; d<Dim; d++)
    res[d] = coord[d]-lhs[d];
  return res;
}

template <class T, int Dim>
inline Vec<T,Dim> operator+(const Vec<T,Dim> &lhs) {
  Vec<T,Dim> res;
  for(int d=0; d<Dim; d++)
    res[d] =coord[d]+lhs[d];
  return res;
}

template <class T, int Dim>
inline Vec<T,Dim> operator*(const Vec<T,Dim> &lhs) {
  Vec<T,Dim> res;
  for(int d=0; d<Dim; d++)
    res[d] =coord[d]*lhs[d];
  return res;
}

template <class T, int Dim>
inline Vec<T,Dim> operator*(double lhs) {
  Vec<T,Dim> res;
  for(int d=0; d<Dim; d++)
    res[d] =coord[d]*lhs;
  return res;
}

template <class T, int Dim>
friend inline Vec<T,Dim> operator*(double lhs,const Vec<T,Dim> &lhss) {
  Vec<T,Dim> res;
  for(int d=0; d<Dim; d++)
    res[d] =lhss[d]*lhs;
  return res;
}

template <class T, int Dim>
friend inline Vec<T,Dim> operator*(const Vec<T,Dim> &lhss,double lhs) {
  Vec<T,Dim> res;
  for(int d=0; d<Dim; d++)
    res[d] =lhss[d]*lhs;
  return res;
}
friend ostream& operator<<(ostream &os, const Vec<T1,Dim> &p) {
    os << "( " << p[0];
    for(int d=1; d<Dim; d++)
      os << "," << p[d];
    os << " )";
    return os;
}
};




inline double norm(double x, int nt)
{
  return std::abs(x);
}

template <class T, int Dim>
inline T norm(const Vec<T,Dim> &lhs, int nt = 2)
{
  if(nt == 2) {
    return sqrt(dot(lhs, lhs));
  } else if(nt == 1) {
    T a = 0;
    for(int i=0; i<Dim; a+=std::abs(lhs[i]), ++i);
    return a;
  } 
  return 0; // never reached
}

template <int Dim>
inline Vec<double,Dim> normalize(const Vec<double,Dim> &lhs) {
  double l = norm(lhs, 2);
  return lhs / l;
}


template <class T, int Dim>
inline T dot(const Vec<T,Dim> &lhs, const Vec<T,Dim> &rhs) {
  T r = 0;
  for(int d=0; d<Dim; d++)
    r += lhs[d] * rhs[d];
  return r;
}



#endif // VEC_H
